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ABSTRACT 

We use a set of numerical TV-body simulations to study the large-scale behavior of the reduced 
bispectrum of dark matter and compare the results with the second-order perturbation theory and 
the halo models for different halo mass functions. We find that the second-order perturbation theory 
(PT2) agrees with the simulations fairly well on large scales of fc < 0.05ft.Mpc~^, but it shows a 
signature of deviation as the scale goes down. Even on the largest scale where the bispectrum can 
be measured reasonably well in our simulations, the inconsistency between PT2 and the simulations 
appears for the colinear triangle shapes. For the halo model, we find that it can only serve as a 
qualitative method to help study the behavior of Q on large scales and also on relatively small scales. 
The failure of second-order perturbation theory will also affect the precise determination of the halo 
models, since they are connected through the 3-halo term in the halo model. The 2-halo term has 
too much contribution on the large scales, which is the main reason for the halo model to overpredict 
the bispectrum on the large scales. Since neither of the models can provide a satisfying description 
for the bispectrum on scales fc ~ 0.1 /iMpc~^ for the requirement of precision cosmology, we release 
the reduced bispectrum of dark matter on a large range of scales for future analytical modeling of the 
bispectrum. 

Subject headings: gravitational Icnsing — dark matter — cosmology: theory — galaxies: formation 



1. INTRODUCTION 

The understanding of the formation, clustering, and 
evolution of galaxies needs to involve the dark matter 
in modern cosmological models. In the current standard 
paradigm, i.e., A cold dark matter (CDM) models, dark 
matter has negligible velocity in the early universe and 
dominates the matter density in the present universe. 
Starting from an initially tiny Gaussian density fluctu- 
ation, the density distribution of CDM would keep its 
Gaussianity throughout the whole evolution in the linear- 
order perturbation (e.g., Peebles 1980; Bernardeau et al. 
2002; Sefusatti & Komatsu 2007). However, the dynam- 
ical evolution of the fluctuations, through the mode cou- 
pling of different scales, would result in non-Gaussianity 
of density distribution on all scales, including the con- 
ventional linear regime. Various statistical tools are 
used in the literature to describe the density distribu- 
tion, and can be further incorporated into the research 
of non-Gaussianity. For example, the two-point corre- 
lation function(2PCF) ^(r) is employed to describe the 
correlation of two objects at a separation r, while the 
three-point correlation function(3PCF) is used to de- 
scribe the next-order statistics in the series of connected 
A^-point correlation functions(NPCFs) (Peebles 1980). 
In the linear-order perturbation theory, a Gaussian field 
is completely determined by the 2PCF. Therefore, the 
3PCF is the lowest-order non-Gaussian quantity in the 
series, which can give us important clues to the proper- 
ties of nonlinear evolution. 

The dark matter 3PCF and its Fourier-space counter- 
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part, the bispectrum, have been studied intensively over 
the years(e.g., Fry 1984; Scoccimarro et al. 1998; Hou 
et al. 2005; Nichol et al. 2006). The theoretical work 
mainly comprises two methods: the nonlinear perturba- 
tion theory(hereafter PT; e.g., Bernardeau et al. 2002, 
and the references therein) and the halo model(hereafter 
HM; Jing & Borner 1998; Ma & Fry 2000; Peacock & 
Smith 2000; Seljak 2000; Scoccimarro et al. 2001; Takada 
& Jain 2003; Yang, Mo, & van den Bosch 2003; Wang et 
al. 2004; Smith et al. 2008). Both methods have been 
applied to and compared with simulations over a range 
of scales(e.g., Scoccimarro et al. 2001; Hou et al. 2005; 
Smith et al. 2008). In addition, as a prior assumption, 
PT is believed to be applicable on large scales where 
the density fluctuation is small enough for a perturba- 
tive method to be useful. It starts to fail as the scale 
goes down to the strongly nonlinear regime. 

The halo model approach, however, is based upon dark 
matter halos. By utilizing halo correlation functions and 
the dark matter distribution in the halos, the halo model 
can then generate the dark matter clustering statistics 
without any consideration of the scales. An advantage 
of the halo model method is that after being combined 
with the models of galaxy distribution in the host halos, 
which is often called halo occupation distribution(HOD; 
e.g., Jing et al. 1998; Berhnd & Weinberg 2002; Yang, 
Mo, & van den Bosch 2003; Zehavi et al. 2004; Zheng 
2004; Zheng et al. 2005), it can produce the overall dis- 
tribution of galaxies, which can be directly compared 
with the observations of large-scale galaxy surveys. So, 
the halo model method can also be used to predict the 
galaxy-biasing relations(e.g., Seljak 2000; Scoccimarro et 
al. 2001; Verde et al. 2002). 

The development of numerical simulations in recent 
years helps test the validity and range of applicability of 
these theories(e.g., Scoccimarro et al. 2001; Fosalba et 
al. 2005; Smith et al. 2008). On large scales, both PT 
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and HM are believed to be in good agreement with the 
simulations. On smaller scales(quasi-linear and nonlinear 
regimes), PT breaks down fairly quickly, while HM may 
still be a good approximation. However, most of the 
previous simulation studies in the literature only measure 
the bispectrum on the scales smaller than that of A; ~ 
0.1 /iMpc^^, due to the limited dynamical range of their 
simulations. 

In this work, we use a set of high-resolution large- 
box cosmological iV-body simulations to accurately de- 
termine the dark matter bispectrum on both largo and 
small scales. The results will be compared with the pre- 
dictions of PT and HM in order to test their validity 
in describing the bispectrum. It is important to deter- 
mine where and how well the theories are valid, since a 
determination of the galaxy bispectrum in observation 
can yield a determination of the galaxy bias factors on 
large scales if the bispectrum of dark matter can be pre- 
dicted (e.g.. Pry 1984; Fry et al. 1993). As we wiU see 
that neither the second-order PT (hereafter PT2; with- 
out loop corrections) nor the HM can provide a satisfying 
description for the dark matter bispectrum on scales of 
k ~ 0.1 hMpc~^ in the standard of precision cosmology, 
we make the bispectrum of dark matter on a large range 
of scales available to interested readers, who want to use 
the data to construct a more accurate model for the bis- 
pectrum. 

This paper is organized as follows. In Section 2, we 
describe some basic formulae and the method used in 
this paper. We introduce our simulations in Section 3. 
We compare the simulations with PT and halo models in 
Section 4.1 and Section 4.2, respectively. We summarize 
our results in Section 5. 

2. METHODS 

In a statistically homogeneous and isotropic universe, 
the density fluctuation of the matter distribution is de- 
fined as (5(x) = (/9(x) — p) /p, with the mean density being 
p = (/9(x)). Due to the convenience in calculating bis- 
pectrum in Fourier space, we transform S{x) to 



Su_ = J rf^xexp(— ik • x)5(x) 



(1) 



By definition (Jk) = 0. So the next two nonvanishing 
terms in the series will be 

(4,<5k.)='5°(ki+k2)P(fci) (2) 
{6i.,S^M = S^iki +k2+ k3)B(fci, fe, fca), (3) 

where (5°(k) is the three-dimensional Dirac delta func- 
tion and P{k) and B{ki, k2, k^) are the Fourier transform 
of the 2PCF and 3PCF, known as the power spectrum 
and bispectrum, respectively. If we extend the pertur- 
bation theory to the second order, it gives (Fry 1984; 
Bernardeau et al. 2002) 



SpT(fcl,fc2,fc3) = i^(ki,k2)PL(fci)PL(fc2) +CI/C, (4) 

ki • k2 

k2 ki 



P(ki,k2) = (l +//) + ■ 



+(1-M) 



kik2 

ki ■ k2 

kik2 



whore Phik^) if^ the linear power spectrum and p = 

3rim^^^^/7 denotes the weak dependence on cosmology. 
It is often convenient to define the reduced bispectrum. 



Q{ki,k2,kz) = 



BpT{ki,k2,k:i) 

PL{kl)PL{k2) + CljC 



(6) 



(5) 



where /ci, ^2, and k^ set the three sides of the triangle 
shape. Two sets of variables are frequently used for Q. 
We define the three variables to be k = ki, u = k2/k\, 
and V = (fca — k2)/ki, where ki ^ k2 ^ k^, u > 1, and 
< w < 1. Therefore, k sets the scale of the triangle, 
while u and v determine its shape. In addition, if we 
change v to a = arccos(ki ■ 'k2)/kik2, (0 < n < n), we 
get the other form as Q{k, u, a), which clearly shows the 
relation between ki and k2. Also, the shape dependence 
is obviously described by u and a. But this parametriza- 
tion has a disadvantage that when considering triangles 
of different scales and different shapes together, we may 
have counted the set of {kx,k2,kz) repeatedly, because 
in this case we only have Q < k^ < k\ + k2- In order to 
avoid repeated counting, we still need to impose a con- 
straint of /c3 > k2 when wo count independent triangles. 
We will use Q{k,u,a') as the preferred parametrization 
when investigating the scale and shape dependence of Q 
in the following sections. 

Although Q(fc, u, a) is generally dependent on the scale 
and shape of a triangle, in the second-order perturba- 
tion theory, if we assume the power spectrum to be a 
power law, Q is then completely independent on the scale 
k. For the equilateral triangle shapes ^1 = ^2 = ^3, 
from Equations (4) and (6) we have Qeq = 1/4 + 3/i/4, 
which is exactly constant. These important properties 
of the reduced bispectrum Q make it a better quan- 
tity than the bispectrum B{ki,k2, ks) to characterize the 
non-Gaussianity of the matter distribution. The shape 
dependence of Q reflects the influence of the gravitational 
instability, so it can be significantly affected by the exis- 
tence of large-scale structiircs such as filaments (Sefusatti 
& Scoccimarro 2005) for small simulation volumes. In 
order to reduce such finite- volume effects, large-box sim- 
ulations are necessary to compute the bispectrum as we 
do in this paper. 

3. AT-BODY SIMULATIONS 

The cosmological model considered here is a canoni- 
cal spatially flat CDM model with the density parame- 
ter = 0.268, the cosmological constant JIa = 0.732, 
the Hubble constant h — 0.71, and the baryon density 
parameter fi^ = 0.045. The primordial density field is as- 
sumed to be Gaussian with a scale-invariant power spec- 
trum cxfc. For the linear spectrum, we generate it from 
the CMBfast code (Seljak & Zaldarriaga 1996), and the 
normalization is set by as = 0.85, where as is the present 
linear rms density fluctuation within a sphere of radius 
S/i-^Mpc. 

We use an upgraded version of the Particle-Particle- 
Particle-Mesh (P^M) code of Jing & Suto (1998, 2002) 
to simulate structure formation in the imiverso. The 
code has now incorporated the multiple level P"^M grav- 
ity solver for high-density regions (Jing & Suto 2000). 
In order to have a large mass resolution range, we run a 
total of 11 simulations with 1024"^ particles in different 
simulation boxes, which we hereafter denote by Z/600, 
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TABLE 1 

Simulation parameters 



boxsize( /i ^Mpc) 


particles 


realizations 


^particle 


600 


10243 


3 


1.5 X lOio?i-iM0 


1200 


10243 


4 


1.2 X lO"?i-iM0 


1800 


10243 


4 


4.0 X 10"/i-iMc7) 



L1200, and L1800 by different box sizes(Table 1) (Jing 
at al. 2007). The simulations were run on an SGI Altix 
350 with 16 processors with OPENMP parallelization in 
Shanghai Astronomical Observatory. 

The different box sizes put different hmits on the de- 
tection scales in the simulations , which is fc > kbox = 
27r/if,orr- For our simulations here, wc have ktox = 
0.0035, 0.0052, 0.01 hUpc^^ for L1800, il200, and L600, 
respectively. Thus, in principle, we expect to deter- 
mine the behavior of the dark matter bispectrum on 
scales larger than k = 0.1 hMpc~^ with good precision. 
The large resolution range used here enables us to check 
the consistency of the results among different simulation 
boxes and to investigate the behavior of dark matter bis- 
pectrum down to small scales. Here, we choose the bin 
scheme as Alogio{k) = 0.1 for fc < 0.1 hMpc^^ and 0.05 
for the larger k in Fourier space. But even for L1800, we 
only have four realizations, which may still be not enough 
to recover the exact variance of the bispectrum on large 
scales. So, we calculate the uncertainty of bispectrum 
under the assumption of a Gaussian field as (e.g., Fry 
et al. 1993; Scoccimarro et al. 1998, 2004; Sefusatti & 
Komatsu 2007) 



1 



Ni 



23 



-Ptot{kl)Ptot{k2)Ptot{k3) 



P,,,{k)=P{k) + —, 



(7) 



where A^i23 is the number of independent triangle con- 
figuration modes in the Fourier space, and Ptotik) is the 
power spectrum that includes the contribution of shot 
noise with Np particles. Compared with Equation (28) 
in Sefusatti & Komatsu (2007), there is no ,si23 factor in 
our Equation (7) since is the number of independent 
triangle modes and the statistically rotational-invariance 
effect has already been taken into account. 

We show the ratios between the variances of the mean 
of the reduced bispectrum averaged over the realizations 
and the corresponding uncertainties of a Gaussian field 
(Equation (7)) in Figure 1, where different panels stand 
for different simulation boxes. The scale k shown in Fig- 
ure 1 is defined to be the longest side of the triangle. 
In addition, to make a complete statistics of the vari- 
ance ratio, we have included all the triangle configura- 
tions in each k bin. On large scales of A; < 0.2 /iMpc~^, 
the Gaussian uncertainties are slightly larger than those 
determined from different simulation realizations, with 
o'sim/o'gaussian ~ 0.9. This is probably due to the fact 
that we have used P{k) determined from the simulations 
to normalize the bispectra (see Equation (6)). As the 
scale goes down to the nonlinear scales where the as- 
sumption of a Gaussian field is not applicable, the ratio 
increases rapidly, as shown in the case of L600. So, for 
our purposes here, the utilization of Gaussian variance 
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Fig. 1. — Ratios between the variances of the mean of the re- 
duced bispectrum averaged over the realizations and the theoretical 
Gaussian uncertainties. Different panels stand for different simu- 
lation boxes. The scale k is defined to be the longest side of a 
triangle. 

on the large scales can, to some extent, compensate for 
still limited number of independent realizations of our 
simulations. 

4. RESULTS 

4.1. Dark Matter Bispectrum 

We mainly present our results with the reduced bispec- 
trum Q{k, u, a). Before we investigate in detail the scale 
and shape dependence of Q, we check the consistency be- 
tween PT2 and the simulations for the equilateral trian- 
gle configurations, which are defined that ki,k2, k^ are in 
the same k-space bins. We show the scale dependence of 
Qeq in Figure 2 where the points in different panels repre- 
sent the different realizations of each simulation set from 
il800 to L600. The solid lines with the error bars are the 
mean values with variances. Here, the variances of the 
mean are obtained from the Gaussian uncertainties in- 
stead of the variances from different realizations in order 
to obtain an accurate error estimation on large scales. Al- 
though they would be underestimated on smaller scales, 
it does not affect our discussion because we mainly focus 
on the large-scale behaviors. The PT2 predictions are 
shown as the dashed lines, with Q^q = l/i+Bn/A = 0.585 
for our cosmological parameters. 



From Figure 2, 



note that different Li, 



simula- 



tions are consistent with each other on nonlinear scales 
of fc > 0.2/iMpc~"^ where Qeq deviates rapidly from the 
PT2 prediction as expected. But the situation is more 
complicated on large scales where large fluctuations from 
realization to realization are displayed, which is also im- 
plied by the Gaussian uncertainties shown in the figure. 
For the case of il800, which clearly shows the trend 
of being constant at fc < 0.1 hMpc^^, the mean of Qeq 
is, however, still about 20% larger than that predicted 
by PT2 at 0.04 /iMpc"^ < fc < 0.1/iMpc"^ We also 
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Fig. 2. — The scale dependence of Qeq- The points in different 
panels represent the different realizations of each simulation set 
from L1800 to L600. The solid lines with the error bars are the 
mean values with variances. The variances here axe obtained from 
the Gaussian uncertainties (Equation?). Also PT2 predictions are 
shown as the dashed lines. 

find such deviations in L1200 at scales of 0.05 hMpc~^ < 
k < 0.1 hMpc~^. For L600 the tendency toward constant 
starts from k ~ 0.06/iMpc~^, but the large fluctuation 
prevents us from telling whether the simulation has a 
similar deviation from the PT2. The jump in the errors 
in L600 at 0.1 /iMpc~^ is caused by our change of Fourier 
space bin scheme at this scale. 

The fluctuation of Q among the different simulation 
realizations on large scales is probably due to the finite- 
volume effect or the scarcity of realizations for each Lbox ■ 
From the fluctuation, we can determine the scales where 
enough independent triangle configurations can be con- 
structed to get an accurate value of Q. As Lbox changes 
from L1800 to L600, the scales where Q can be deter- 
mined are k > 0.03/;,Mpc"\ k > 0.04/iMpc"\ and 
k > 0.08 /iMpc^^, respectively if we require the errors 
are less than about 30%. Thus the finite-volume effect 
is still very important at large scales of < 0.1 hMpc~^ 
even for simulations with a size as large as that of L600. 
It would be extremely cautious to calculate the dark mat- 
ter bispectrum on scales of < 0.1 hMpc~^ when using 
simulations with a size much smaller than that of L600. 
Moreover, to make an accurate measurement on scales 
up to A; ~ 0.01/iMpc~^, we still need to involve more 
realizations or simulations of an even larger box size. 

In Figure 3, we show the scale and shape dependence 
of Q{k, u, a) from large to small scales with u = fc2/fci = 
2. Different points denote different Lbox results, with 
circles for L1800. squares for L1200, and triangles for 
L600. The error bars are the variances of the mean from 
the different realizations. The dashed lines represent the 
PT2 predictions. As we see from the figure, different Lbox 
simulations are consistent with each other over all the 
scales. So, different simulation resolutions do not affect 



tlie measurement of Q over all the scales here. The PT2 
predictions only agree well with the simulations on scales 
of fc < 0.05/iMpc~^ (where the measurement errors are 
also large). Starting from k ~ 0.06/iMpc"\ PT2 has 
already deviated from the simulation predictions, which 
indicates that the nonlinear effect may be important even 
on these large scales. 

We check the influence of nonlinear corrections by sim- 
ply replacing the linear power spectrum Phik) in Equa- 
tion (4) and Equation (6) with the nonlinear power spec- 
trum PNL{k) derived directly from the simulations. The 
reason to do this is that in real observations, both the bis- 
pectrum and the power spectrum are nonlinearly evolved 
quantities. The results are shown as the solid lines in 
Figure 3. This nonlinear correction can improve some- 
what the match between the simulations and the PT2 at 
scales between 0.05/iMpc~^ and O.l/iMpc"^, but does 
not remove the mismatch completely. 

Despite the overall agreement of PT2 with the sim- 
ulations on large scales for many different a configura- 
tions, PT2 still overestimates the true Q{k,u,a) for co- 
liner triangle configurations, as shown clearly in the top 
left panel of Figure 3. This phenomenon was also discov- 
ered by many other authors(e.g., Scoccimarro et al. 1999; 
Smith et al. 2008). It may be explained by the fact that 
these large-scale structures tend to be filamentary rather 
than spherical in the simulations. For configurations near 
that of the isosceles triangles where ks = k2 = 2k\ and 
a/-K = 0.58, there is, however, a trend that PT2 un- 
derestimates the bispectrum from k = 0.06 /iMpc""'^ to 
0.1/iMpc~^. Such underestimates are not due to the 
shot-noise effect in the simulations. We have corrected 
for the shot-noise effect when obtaining Q. Because we 
have 1024^ dark matter particles for each simulation, the 
shot noise is negligible for the final results. So, such an 
inconsistency demonstrates the defect of PT2 on these 
large scales. In the framework of perturbation theory, 
to describe the simulations better, we might have to in- 
corporate higher-order corrections, such as the one-loop 
correction(Scoccimarro & Frieman 1996; Bernardeau et 
al. 2002). Involving such corrections may change the 
shape of Q even on large scales, i.e., to make Q smaller 
for colinear triangles and larger near isosceles triangles, 
consistent with our results (e.g., Bernardeau et al. 2008). 

4.2. The Halo Model 
4.2.1. Formalism 

To better understand the nonlinear effects in the re- 
duced bispectrum, we use the halo model method to the- 
oretically study the properties of Q{ Jing et al. 1998; Jing 
& Borner 1998; Ma & Fry 2000; Peacock & Smith 2000; 
Seljak 2000; Scoccimarro et al. 2001; Berhnd & Wein- 
berg 2002; Cooray & Sheth 2002; Takada & Jain 2003; 
Wang et al. 2004; Smith et al. 2008). The essence of 
the halo model is to construct the distribution of parti- 
cles(dark matter or galaxies) by fixing the distribution 
of the clumps (dark matter halos) and the distribution of 
the particles within the clumps on the assumption that 
all the particles reside in the clumps. So, it is a halo- 
based method and the changes in the properties of halos 
will affect the final distribution of the particles. 

Such a hierarchical method is very easy to apply given 
the information of the different model components. The 
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Fig. 3. — Scale and shape dependence of Q{k, u, a) from large to small scales with u s 2, i.e., k2 = 2ki. Different points denote the results 
for different Lf,g^, with circles for L1800, squares for L1200, and triangles for L600. The error bars are the variances on the average of the 
different realizations. The dashed and solid lines represent the PT2 predictions with the linear power spectrum Phik) and the nonlinear 
power spectrum PjsrLik), respectively. 



dark matter halo distribution is generally constrained by 
the halo mass function n{m), which describes the dark 
matter halo number density n as a function of halo mass 
TO. The particle distribution within the halos is defined 
by the dark matter halo density profile p{r\m). So, the 
particle particle correlations can be obtained from the 
halo-halo correlations. By assuming a biased distribu- 
tion of halos relative to the underlying mass, we can 
describe the halo halo correlations with the halo bias 
parameters h{m){i = 1,2, . . .)(Mo, Jing & White 1997; 
Sheth & Tormen 1999). 

Following the notation described by Cooray & Sheth 
(2002), the dark matter power spectrum can be decom- 
posed into the 1-halo and 2-halo terms as 

P{k) = Pihik) + P2h{k). (8) 
These terms are 

Pih{k) = Mo2{k,k) 

P2h{k) = [Mu{k)fPL{k), (9) 

where 



Mij{ki, . . . ,kj)= J rfmn(TO) J bi{m) 

x[u{ki\m) . . .u{kj\m)]. (10) 



and &o = 1- 

Here, p is the mean density of the universe and u{k\m) 
is the normalized Fourier transform of the dark matter 
halo density profile p{r\m) truncated at the virial radius 

/7 1 N r^<^'~ osinfcr p(r\m) ,,,, 

u(k\m) = / dr Awr^ — — ^. (11) 

Jo kr m 

As A; — *■ 0, u{k\7n) — > 1 and by definition, we have 

/Tfl 
dmn{m)—bi{m) = 1, (12) 

and so on large scales P-inik) Phik), as expected. 

The bispectrum is similarly decomposed into the con- 
tributions from the 1-halo, 2-halo, and 3-halo terms, 

Bi23 = Bih + B2h + B^h, (13) 

where 

Bih = Mo^{ki,k2,kz), 

B2h = Mn{ki)Mi2{k2,k3)PL{ki)+cyc., 

Bsh = Mil (fci)Mii (/c2)Mii {k^)BpT 

+ [Mn{ki)Mn{k2)M2i{k3)PL{ki)PL{k2) + eye] . 

(14) 
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Again we have B^h Bpx as fc ^ 0, which means that 
the 3-halo term will return to that of the PT2 predictions 
on large scales. 

The reduced bispectrum Q is then defined as 



B 



123 



P{kr)P{k2) + eye 



(15) 



Note the difference from Equation (6) in the denominator 
where the linear power spectrum Pl {k) is replaced by the 
nonlinear power P{k) determined in the halo model. 

In our model, dark matter halos are defined as objects 
with a mean density A^j^ times that of the background 
universe(Bryan & Norman 1998), where A„ir ~ 361 
for our cosmology parameters, and their density distri- 
bution follow the Navarro Frcnk White (NFW) profile 
(Navarro, Frenk, & White 1996). The concentration pa- 
rameter c(rn) in the density profile is given by the re- 
lation c(m) = co{m/m^)^ , where cq = 9, /? = —0.13, 
and TO* = 4.8 x 10^^ H^^Mq is the nonlinear mass 
scale(Bullock et al. 2001). We also generate the linear 
power spectrum from CMBfast here. 

For the halo mass function(MF), first we consider 
the two commonly used analytical models of Press & 
Schechter (1974)(PS) and Sheth & Tormen (1999)(ST). 
They are defined by 



n(m)dm = —/Mdv, 

TO 



(16) 



where v = 6c/(j{m), with 5c = 1.69, the spherical 
collapse threshold of the linear overdensity 5. a{m) 
is the linear rms fluctuation within spheres of radius 
R = (3TO/47rp)i/3, 



I 



dfc k^PL{k) 



W'^{kR), 



(17) 



k 27r2 

where W{x) = 3(sinx — xcosx)/x^ is the Fourier trans- 
form of the top- hat window function and W{x) 1 as 
X ^ 0. 

The function /{ly) for PS and ST can be generalized 
into the following form; 




13 14 
Log,„{m) [ Ujh ] 

Fig. 4. — Three mass functions compared with our simulations. 
Different points stand for different I/joa, simulations and different 
lines for different mass functions 

/■^^ TO 

1 — / dmn{m)—bi{m) 



(20) 

The upper mass limit M2 is usually set by the maximum 
halo mass in the simulations. Here, we just use the max- 
imum halo mass seen in L1800 as M2 = 6 x 10^^ /i~^Mq. 

For the halo bias parameters, we use the results of Mo, 
Jing & White (1997) and Sheth & Tormen (1999). We 
only need to consider the first two halo bias parameters: 



where 



hi{m) = l + ex + Ei, (21) 
h2{m) = 2{l + a2){ei+Ei) + e2 + E2, (22) 



./(.) = 2A^ [I + {a.yP] exp (^-^) (18) 

where A = 1,0.322, p = 0,0.3, and a = 1.0.707 for PS 
and ST, respectively. 

Jenkins et al. (2001) have also derived the halo mass 
functions for different halo definitions of FOF halos and 

spherical ovcrdensity(SO) halos. Here, we only use their 
FOF halo mass function to match our halo definition 

i//(i/) = 0.307exp(-| ln(i//5c) -|- OMf-^"^), (19) 

with —0.9 < \n{i//5c) < 1-0, corresponding to about 
4 X IO^/i-^Mq < to < 4 X IO^^/i-^Mq in our simu- 
lations. Although this mass function is derived within 
a finite halo mass range and does not satisfy the den- 
sity normalization requirement, we can still impose the 
normalization constraint of Equation (12) by only inte- 
grating over a limited mass range, i.e., 

f°° TO Z"*^^ TO 

Jo P J Ml P 



av — 1 av 5, 

ei = 7 , £2 = - 3), 

<5c 6^c 

2p/5c E2 l + 2p 
1 -I- {av^y El 5c 

and 02 = -17/21. (23) 

As for the Jenkins et al. (2001) MF(hereafter JMF), we 

also use the bias parameters corresponding to the PS 
MF. Changing to the other halo bias model does not sig- 
nificantly affect the predictions for the power spectrum 
or bispectrum. 

We compare the three mass functions with our simu- 
lations for 'w?"n['m) / p in Figure 4, where different points 
stand for different Lbox simulations and different lines 
for different mass functions. The error bars arc deter- 
mined from the scatter among the different realizations. 
The fiuctuations shown in the low mass ends of differ- 
ent Lbox simulations are caused by our lowest mass halo 
definitions that we set 10 particles as the minimum re- 
quirement for a bound halo. Except for the PS MF, 
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Fig. 5. — Ratio between the nonlinear power spectrum P{k) and 
the hnear power Pj. (k) for the halo model predictions with different 
MFs, compared with those of the simulations. 

the other MFs have similar halo number densities at the 

high mass end, though JMF is somewhat smaller than 
ST MF for TO > 10^^ H-^Mq. The figure indicates that 
PS MF has much fewer large halos while retaining too 
many small mass halos. The ST MF is consistent with 
the simulations for the high mass and low mass halos, 
but it is substantially lower in the range of lO^'^ H^^Mq- 
lO^^ft-~^M0, which is quite important for the scales we 
consider here. JMF seems to agree with the simulations 
better for these intermediate masses, though in a detailed 
comparison we find that it is still lower than the simula- 
tions for TO > 10^^ h~^MQ. 

In our simulations, about 20%- 30% of the total halo 
particles reside in halos with m > 10^^ /j^^Mq and about 
53%-80% are in halos of m > lO^^ H-^Mq. The fluctua- 
tions are due to the different resolutions of different Ltox 
simulations. Therefore, the differences among different 
MFs in the mass range of 10^^ /i^^Mq-IO^^ H^'^Mq wiU 
strongly affect the final results at least in the intermedi- 
ate scales as shown in Figure 5. Although we do not have 
a mass function exactly consistent with the simulations, 
we can still compare the results of Q for these three MFs 
to study the effect of different halo mass distribution. 

4.2.2. The Halo Model Results 

After specifying the MF and the halo bias, we can 
theoretically determine the power spectrum and the bis- 
pectrum of the dark matter. In Figure 5, we show the 
ratio between the nonlinear power spectrum P{k) and 
the linear power PL{k) for halo models of different MFs, 
compared with those of the simulations. The wiggles on 
the large scales are caused by the fact that the Bary- 
onic Acoustic Oscillations (B AO) features are suppressed 
in the nonlinear evolution relative to the linear evolution. 
The nonlinear power spectrum from the simulations is a 
bit smaller than the linear power on large scales (corre- 
sponding to the first BAO peak), and it rapidly increases 



when k > 0.1/iMpc^ . The halo models for different 
MFs converge on large scales where the 2-halo term dom- 
inates. But they are still larger than the linear power 
spectrum, because the 1-halo term still has a nonvanish- 
ing contribution to P{k) on these large scales. Compared 
with the simulations, the power spectra of the halo mod- 
els arc larger than those of the simulations, and the wig- 
gles in A: > 0.1 hMpc~^ where the 1-halo term becomes 
important are not as evident as in the simulations. For 
scales of fc > 0.2 hMpc~^ , all the halo model predictions 
are substantially lower than those of the simulations. On 
such scales, the large mass halos play an important role 
in the 1-halo term. And yet we have not corrected for the 
halo boundary and exclusion effects(e.g., Takada & Jain 
2003) in the halo models. The effect of choosing an MF 
is relatively small, as the predicted power spectra based 
on the three MFs are very similar. The simple analytical 
prescription of the halo model, instead of its ingredients 
(the MF, the halo density profile, and the halo bias fac- 
tors) , may be the main cause leading to the discrepancies 
between the models and the simulations shown in Figure 
5. 

In Figure 6, we show the halo model predictions for 
the reduced bispectrum Q{k,u,a). As in Figure 3, the 
points are for the simulations. The dotted, dashed, and 
solid lines show the halo model predictions for PS, ST, 
JMFs respectively. Compared with the PT2 predictions 
in Figure 3, the halo models are a bit larger on large 
scales, as shown in the top panels of Figure 6. Differ- 
ent halo model predictions converge on the large scales 
where the 3-haIo term dominates. But the predictions of 
the JMF model are somewhat larger than the others. It 
is generated by the fact that the 2-halo term contribu- 
tion is much higher for JMF. Such differences are shown 
clearly in the intermediate scales of 0.06 0.1 ftMpc""'^ in 
the middle panels of Figure 6. On these scales, the halo 
models capture the shape of reduced bispectrum in the 
simulations rather well compared with the PT2 predic- 
tions, though the amplitudes are still higher than ex- 
pected. The result for the PS MF seems to fit with the 
simulation bispectrum better near the isosceles triangle 
shapes. For the quasi-linear and nonlinear scales of 0.2- 
0.5 hMpc~^ in the bottom panels of Figure 6, the differ- 
ences between the MFs show up distinctly. The model 
predictions for the PS MF are in much better agree- 
ment with the simulations. Even at the highly nonlinear 
scale of A; ^ 0.5/iMpc^^, it can provide a fairly good 
estimate of Q. And as k increases, the reduced bispec- 
tra of the simulations approach a constant more slowly 
than the halo models predict. Moreover, the turnover 
feature in the shape dependence, appeared in the halo 
model predictions for the ST and JMFs at the scale of 
k ~ 0.5/iMpc^^, is not seen in the simulations. 

The better agreement on the reduced bispectrum be- 
tween the halo model predictions for the PS MF and the 
simulations on scales of A: > 0.2/iMpc~^ does not mean 
that PS MF is a better description for the mass function 
in halo models. Since the reduced bispectrum is the ratio 
of bispectrum to a sum of power spectrum products, the 
better agreement achieved with the PS MF is largely a 
coincidence, because the power spectra are fitted worse 
on relevant scales by the PS model(as shown in our Fig- 
ure 5). In fact, the ST MF and JMF better describe 



8 




Fig. 6. — Same as Figure 3, but for the halo model predictions with different MFs. The points are for the simulation results. The dotted, 
dashed, and solid lines show the halo model predictions with PS, ST, JMFs respectively. 



the mass function of halos in our simulations (Figure 4). 

We also note that our result is consistent with Fosalba 
et al. (2005), who argued that a better agreement in the 
reduced bispectrum between the halo model and simula- 
tions can be achieved by imposing a mass cutoff for halos 
at the high mass end. But, we do not think that impos- 
ing the mass cutoff is physical in our study, since our 
simulation box is big enough that the halos at high mass 
end of ST MP (or JMF) are fairly well sampled in our 
simulations. Therefore, the failure of the halo model for 
the ST MF or JMF imphes the inherent difficulty of the 
halo models to accurately predict the power spectrum 
and the bispectrum. 

To investigate the c;fFc!cts of choosing different MFs 
and the contributions of different halo terms, we show 
in Figure 7 the scale and shape dependence of Bih / B123 
of each halo term by changing k and a while fixing 
H = ^2/^1 = 2. We show from top to bottom for 
three triangle shapes including the colinear and isosceles 
triangles, specified by different a values. At the scale 
k = 0.01 /iMpc~^, in each panel from top to bottom are 
the 3-halo, 2-halo, and 1-halo term contributions. Dif- 
ferent lines denote the halo model predictions for dif- 
ferent MFs as in Figures 5 and 6. On large scales of 
k < 0.1/iMpc~^, the 3-halo term dominates but 2-halo 
term still has a contribution of about 5%-20% in these 
halo models. Comparing Figure 3 and Figure 6, it seems 
that the 2-halo term contributions in the halo models 



are higher than simulations predict. Furthermore, we 

note that in the determination of the reduced bispec- 
trum Q{k,u,a) of the halo models, the linear power 
spectrum Pbik) in the denominator of Equation (6) in 
PT2 is replaced by the nonlinear power P{k) which is 
larger than Phik) over all the scales shown in Figure 
5. Therefore, the 2-halo term is actually necessary to 
fit with the simulations, otherwise only using the 3-halo 
term(_B3;i BpT on large scales) will make Q{k, u, a) of 
the halo models much smaller than those of simulations. 
For the different halo term contributions, the differences 
among choosing different MFs are very small on large 
scales. The; PS and ST MFs have nearly the same halo- 
term contributions at these scales. This indicates that 
the bispectrum on these scales is not sensitive to the 
choice of MFs. 

When k > 0.1 hMpc~^, we see from Figure 7 that the 
2-halo term is much more important for the isosceles tri- 
angle shapes than for the colinear ones. It has almost the 
same contribution as the 3-halo term at A; ~ 0.1 hMpc~^ 
for isosceles shapes, while the 3-halo term still dominates 
in the colinear configurations. Thus, it explains the fail- 
ure of PT2 which only includes the 3-halo term for the 
isosceles triangles, shown in the middle panels of Fig- 
ure 3. On these quasi-linear scales, the 1-halo term also 
starts to become important and its contribution differs 
for different triangle shapes, which clearly shows the in- 
fluence on the shape dependence of Q{k,u,a). Again, 
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Fig. 7. — Scale and shape dependence of each halo term con- 
tribution Bj/,/Bi23 with /i = fe2/fei = 2, as well as the effects of 
the different MFs. We show from top to bottom for three trian- 
gle shapes including the colinear and isosceles triangles, specified 
by different a values. At scale k = 0.01 /iMpc~^, in each panel 
from top to bottom are the 3-halo, 2-halo, and 1-halo term contri- 
butions. Different lines denote the results for different MFs as in 
Figure 5 and Figure 6. 

the 1-halo term is most important for the isosceles trian- 
gle configurations. We can therefore infer that the non- 
linear effect is more significant for such triangle shapes. 
Meanwhile, the differences also emerge on these quasi- 
linear scales when different MFs are used. We note from 
Figure 6 that the reduced bispectrum for the ST MF 
case increases rapidly in the quasi-linear and nonlinear 
regime. From Figure 7, we can see that it is caused by 
the rapid increase in the 1-halo term contribution. Even 
in the nonlinear scale of fe = 0.5/iMpc^^, where the 1- 
halo term is dominant, the 2-halo contribution is still not 
negligible and the 3-halo term also has a nonvanishing 
contribution for some triangle configurations. 

Despite the fact that the halo model can provide a 
qualitative explanation for Q{k, u, a) down to the quasi- 
linear scales, there still exist unavoidable problems with 
the halo model if we consider the possibility that the 
halo model can match the power spectrum P{k) of the 
simulations extremely well. We take k = 0.1 hMpc~^ 
for example. We see from Figure 5 that the nonlinear 
power P{k) of the simulations is almost equal to the lin- 
ear power Phik), only with a small deviation at this scale. 
Then wc can replace P{k) in Equation (15) with Pilk). 
Thus, the reduced bispectrum Q{k,u,a) must be larger 
than QpT in PT2, because the 3-halo term has already 
included Bpt (at this scale Mii(fc) is still very close to 
1) but the 2-halo term can even add an additional sizable 
contribution. For the case of colinear configurations with 
a — > 0, in Figure 3, the PT2 prediction is much larger 
than those of the simulations. Therefore, the halo model 



prediction with the 2-halo term included would give an 
even larger bispectrum. This apparent contradiction in 
the predictions for P{k) and Q{k, u, a) between the halo 
models and the simulations actually comes from the pre- 
scription of the halo model. Because in the halo model, 
it includes the PT2 predictions in the 3-halo term in 
Equation (14). The failure of PT2 will then lead to the 
failure of the halo model if its predictions are larger than 
the simulations. Moreover, the 2-halo term and even the 
1-halo term are not negligible at relevant large scales. 

5. DISCUSSIONS AND CONCLUSIONS 

We use a set of numerical A^-body simulations to 
study the large-scale behavior of the reduced bispectrum 
Q{k, u, a) and compare it with the analytical predictions 
from both the second-order perturbation theory (PT2) 
and the halo models. We investigate the scale and shape 
dependence of Q and find that although PT2 agrees 
fairly well with the simulations on the large scales of 
k < 0.05 /iMpc^^, it still shows a signature of deviation 
as the scale goes down. Even on the largest scale where 
the bispectrum can be explored with our simulations, 
the inconsistency between PT2 and the simulations ap- 
pears for the colinear triangle shapes. The better agree- 
ment of PT2 with the simulations after the nonlinear 
power spectrum is used for the PT2 implies the impor- 
tance of nonlinear effects even on these large scales of 
k < 0.1 /iMpc~^, which is a supplement to the previous 
studies(e.g., Bernardeau et al. 2002; Smith et al. 2008). 
Then one of our main conclusions is that higher-order 
corrections (e.g., loop corrections) are necessary on large 
scales of A: = 0.05-0.1 hMpc~^ for the perturbation the- 
ories to predict Q at an accuracy of percent level. 

An alternative method of predicting the reduced bis- 
pectrum is through the halo model with several basic 
ingredients. After a detailed comparison between the 
predictions and our simulation results, however, we are 
disappointed to find that halo models can at best serve 
as a qualitative method to help study the behavior of Q 
on large scales and also on relatively small scales. But 
because of the nonnegligible two-halo term contribution, 
the halo model predictions of the reduced bispectrum 
are always higher than the results of the simulations. 
No halo model can be brought into agreement with the 
bispectrum on nonlinear scales unless a reasonable mass 
function is adopted. 

We have carefully and accurately measured the bis- 
pectrum for a large range of scales in the sim- 
ulations. The results could be useful for fu- 
ture analytical modeling of Q. Therefore, a ta- 
ble of Q from our measurement is made available at 
http://www.shao.ac.cn/mppg/guo/paper/dmbisp.html. 
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